function [Bstd Coor]=Mercury_cusp(IMFxyzT)
 
% -------------------------------------------------------------------------------------------------------------------------------------------------------- 
% An empirical model of Mercury's magnetospheric cusp activity 
% -------------------------------------------------------------------------------------------------------------------------------------------------------- 
% Version 1.0   
% An cusp indicator as a function of solar wind variables and Mercury solar orbital phase is derived from MESSENGER magnetometer data and proton flux data collected between 23 March 2011 and 17 March 2014.
% Details are documented in http://doi.wiley.com/10.1002/2016JA023687
% --------------------- 
% Reference
% --------------------- 
% He, M., J. Vogt, D. Heyner, and J. Zhong (2017), Solar wind controls on Mercury's magnetospheric cusp, J. Geophys. Res. Space Physics, 122, http://doi.wiley.com/10.1002/2016JA023687
% He M., Vogt J. (2018) Empirical Modeling of Planetary Magnetospheres in Response to Solar Wind Dynamics Using EOF Analysis and Multivariate Linear Regression. In: L��?hr H., Wicht J., Gilder S.A., Holschneider M. (eds) Magnetic Fields in the Solar System. Astrophysics and Space Science Library, vol 448. Springer, Cham, http://link.springer.com/10.1007/978-3-319-64292-5_7
%------------------------
%Contents
%------------------------
% Mercury_cusp.m 	:  the model code  v1.0
% CuspDistCoef.mat	: coefficients for the model contained in one variable CuspDist:
% CuspDist.coor.lat and CuspDist.coor.LT, each of which sizes 250x1, define latitude and local time of  250 discrete points in Mercury solar orbital (MSO) coordinates system. 
%  The other coordinates are define as,
% >> CuspDist.coor.xy_comp =cosd(CuspDist.coor.lat).*exp(1i*mod((CuspDist.coor.LT+12)/24,1)*360/180*pi);
% >> CuspDist.coor.y=real(CuspDist.coor.xy_comp );
% >> CuspDist.coor.x=imag(CuspDist.coor.xy_comp );
% CuspDist.Coef_1BxByBzObitObbxObbyObbz sizes 8x250, whoes first dimension  corresponds to are the eight coefficients (two scales pluse two 3d vectors) defined in Equation (2) in http://doi.wiley.com/10.1002/2016JA023687, in the oder:  
% in oder to,
% 1) \beta_0
% 2) \beta_IMFx
% 3) \beta_IMFy
% 4) \beta_IMFz
% 5) \tilde{\beta_t0}
% 6) \tilde{\beta_t0}_IMFx
% 7) \tilde{\beta_t0}_IMFy
% 8) \tilde{\beta_t0}_IMFz 
%------------------------
%Communication to:  
%------------------------
% Maosheng He
% Schlossstrasse 6
% 18225, Kuehlungsborn, Germany
% +49 (0) 38293 68 240
%  hmq512@gmail.com /he@iap-kborn.de
% July, 01, 2019
%%  
%-----------------------
% Syntax
%------------------------
% [Bstd Coor]=Mercury_cusp(IMFxyzT) returns the cusp activity index Bstd on the grid Coor
%------------------------
% Inputs
%------------------------
% IMFxyzT: sizes nx4, whos 2nd dimension corresponds to IMF Bx, By, Bz components and  the number of days since January 0, 0000
%------------------------
% Outputs
%------------------------
% Bstd: activity index (nT);sizes nx250,
% Coor: the coordinates of the 250 points
%------------------------
%Examples
%------------------------
% [x,y,z,t]=ndgrid(-20:5:20,-20:5:20,-20:5:20,0:8:85);
% IMFxyzT=[x(:),y(:),z(:),t(:)];
% [Bstd Coor]=Mercury_cusp(IMFxyzT);
% scatter(Coor.x,Coor.y,500,Bstd(1,:),'.');axis equal;axis ij
%% 
% Copyright (c) 2017, Jacobs University Bremen 
% All rights reserved. 
% Redistribution and use in source and binary forms, with or without modification, are permitted provided that the 
% following conditions are met: 
%  1. Redistributions of source code must retain the above copyright notice, this list of conditions and the 
% following disclaimer. 
%  2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the 
% following disclaimer in the documentation and/or other materials provided with the distribution. 
%  
% THIS SOFTWARE IS PROVIDED BY JACOBS UNIVERSITY BREMEN "AS IS" AND ANY EXPRESS 
% OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF 
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO 
% EVENT SHALL JACOBS UNIVERSITY BREMEN BE LIABLE FOR ANY DIRECT, INDIRECT, 
% INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
% LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE 
% OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF 
% ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 

Obit_Phasor=exp(IMFxyzT(:,end)/87.969*2*pi*1i);
XX(:,1)=Obit_Phasor*0+1;
XX(:,2:4)=IMFxyzT(:,1:3);clear IMFxyzT
XX(:,5:8)=XX(:,1:4).*Obit_Phasor(:,ones(1,4));
load CuspDistCoef.mat
Bstd=XX*CuspDist.Coef_1BxByBzObitObbxObbyObbz;
Coor=CuspDist.coor;


